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Free-floating planets (FFPs)|are planetary-mass objects that are not bound to host stars. First 


discovered in the 1990s, their nature and origin are still largely unconstrained because of a 
lack of large homogeneous samples enabling a statistical analysis of their properties. To date, 
most|FFPs| have been discovered using indirect methods; micro-lensing surveys have proven 
particularly successful to detect these objects down to a few Earth masseg? However, the 
ephemeral nature of micro-lensing events prevents any follow-up observations and individual 
characterisation. Several studies have identified in young stellar clusters4 and the 
Galactic field?! but their samples are small or heterogeneous in age and origin. Here we 
report the discovery of between 70 and 170 (depending on the assumed age) in the 
region encompassing and the closest young OB 
association to the Sun. It is the largest homogeneous sample of nearly coeval|FFPs|discovered 
to date. We found an excess of by a factor of up to seven compared to core-collapse 
models predictions®®, demonstrating that other formation mechanisms may be at work. We 
estimate that ejection from planetary systems might have a contribution comparable to that 
of core-collapse in the formation of Therefore, ejections due to dynamical instabilities 


in giant exoplanet systems must be frequent within the first 10 Myr of a system’s life. 


To date, most exoplanets have been detected through radial velocity or photometric modu- 
lations induced in their host starsW2l As such, the vast majority of known exoplanets are gravita- 
tionally bound to stars. However, several [free-floating planets] have been discovered over 
the last two decades in astro-photometric surveys of nearby star-forming regions ELAES, young 
association? the solar neighbourhood2P4 and in gravitational micro-lensing surveys of the 


Galactic field These ultra-faint objects are incapable of sustaining nuclear fusion and steadily 


fade in time, making them easier to observe when they are very young. [FFPs]are compact objects 
of less than about 13 Jupiter masses that are not bound to a star or brown dwarf but rather wander 
among them. At present, four scenarios have been proposed for the formation of these extreme ob- 
jects: a) a scaled-down version of star formation via core-collapsé25!2% b) within a protoplanetary 
disc, either like gas-giant planets through core accretion or like companions through gravita- 
tional fragmentation of massive extended discs2129 followed by ejection by dynamical scattering 
between planets in both cases=% c) as aborted stellar embryos ejected from a stellar nursery before 
the hydrostatic cores could build up enough mass to become a stat] and d) through the photo- 
erosion of a prestellar core by stellar winds from a nearby OB sta?! While direct observational 
evidence confirms that these different processes are all at work#314] we still do not understand their 


relative contributions to the overall [FFP] population. 


Here we present a search for|FFPs|in the 171 deg? region occupied mainly by [Upper Scorpius 
(USC) and [Ophiuchus| (Oph). We selected an elliptical area centred in (RA = 243.5°, Dec = 


—23.1°) with a semi-major axis of 8.5° in RA and a semi-minor axis of 6.4° in Dec where the spatial 
and temporal coverage of the observations is the best. This large complex represents the perfect 
hunting ground to search for young and nearly coeval [FFPs] thanks to its proximity (120-145 pc) 
and youth (1-10 My BES), We combined our ground-based observations in the optical and infrared 
with wide-field images available in various public archives (see methods). We processed and 
analysed a total of 80 818 individual wide-field images acquired with 18 different cameras over the 


past 20 years to obtain a final catalogue, the|Dynamical Analysis of Nearby ClustErs (DANCe)2 


catalogue, containing positions, proper motions and multi-wavelength photometry (grizyJH K) 


for more than 26 million objects. We complemented the DANCe catalogue with the astrometry 
and photometry of the and Hipparcos! catalogues. We used 
this dataset to compute membership probabilities to[USC]and [OphJusing a probabilistic model of 
the distribution of the observable quantities in both the cluster and background populations (see 
methods). To identify the [FFPs] we used the parameter space that contains the largest amount of 


parameters and sources with complete information, namely proper motions and ¿J H K photometry. 


We identified 3455 high probability candidate members in the area covered by our study, 
including between 70 and 170[FFPs|depending on the age assumed for the region. Approximately 
20% of the members are new compared to previous studies "4240 and this proportion increases 
to 75% in the planetary mass domain. It is the largest and most homogeneous sample of 
identified using direct images so far, and it constitutes an excellent benchmark to test star and 
planet formation theories. Figure [I]shows the area covered by our survey in the optical and radio 
wavelengths. The members are distributed all over the region. The majority of members (and[FFPsp 


lay in the area of[USC| where the extinction by the interstellar medium is lower (Ay < 1 mag). 


We used our comprehensive membership analysis to study the origin of the population 
in[USC]and[Oph| We computed the number density of members as a function of their masses. This 
“mass function” of the region constitutes a fundamental constraint for formation theories because 
different mechanisms predict different relative abundances of stellar, sub-stellar and planetary- 
mass objects. The transformation from observed luminosities to masses requires knowing the age 


and using evolutionary models. The age spread of the regiones as well as the complex and 


overlapping spatial and kinematic distributions of the different coexisting populations, prevented 
us from disentangling the various groups and assigning ages to individual objects. Instead, we 
chose to assign ages of 3, 5 and 10 Myr for the entire sample and assumed that the underlying real 
mass distribution must be included between these borderline cases. We used the entire spectral 
energy distribution available for each object and the theoretical evolutionary models“! to infer 
the individual mass of each member (see methods). Figure[2|shows colour-magnitude diagrams of 
the members in two different spaces. We overplotted the theoretical isochrones at the extreme ages 
mentioned above to illustrate the uncertainties in the mass inherited from the lack of individual 


precise ages. 


The J apparent magnitude distribution (Fig. [3| upper panel) is a direct product of the ob- 
servations, and thus, it is not affected by the uncertainties and errors in the transformation from 
luminosities to masses. Thanks to improved statistics achieved with our larger sample we unam- 
biguously identified a dip at planetary masses (apparent magnitude J ~ 17.8 mag, corresponding 
to absolute magnitude My ~ 12 mag and masses 7-13 Min depending on the age assumed) 
which was also reported in volume-limited samples in the solar neighbourhood“! Are we seeing 
changes relative to the formation mechanisms? Could this be the “real” frontier between brown 
dwarfs and {FFPs) To our knowledge, the origin of this dip remains unknown and more studies are 
needed to answer these questions. However, the presence of this feature at the young ages of 
and [Oph] (1-10 Myr) suggests that it must be the result of the formation and/or early evolution of 


these objects. 


Figure [3| (middle panel) shows the mass function at 5 Myr (which can be considered an in- 
termediate age of the entire complex); the mass functions at 3 and 10 Myr define the upper and 
lower confidence intervals. A number of details and features clearly significant in the magnitude 
distribution (in particular the dip mentioned above) are blurred and lost in the mass distribution 
likely due to the uncertainties related to the transformation from luminosity to mass. Our magni- 
tude and mass distributions show a rich population of We evaluated the contamination rate 
in this sample, mostly due to background reddened giant stars and background Galaxies, using two 
different methods. First, using synthetic data! we estimated a contamination rate in the planetary 
mass domain smaller than 4%. Second, using sources for which both DANCe and|Gaia DR2|mea- 
surements are available (in the range 9 < J < 14 mag, hence above the planetary mass regime) 
and assuming that the [Gaia DR2] sample represents the ground truth, we estimated a contamina- 
tion rate of approximately 8%. The real contamination rate of [FFPs]is likely in between these two 
values and should be confirmed by follow-up spectroscopic observations. The number of 
reported in our analysis must constitute a lower limit of the actual total number of since our 
analysis is expected to miss the objects most affected by extinction (Ay = 3 mag), as well as 
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objects displaying a large near-infrared excess related to the presence of circumstellar material. 


The fraction of [FFPs| in our sample, meaning the relative proportion of to stars and 
brown dwarfs, is 0.04519:953, where the uncertainties come from the uncertainty on individual ages 
(between 3 and 10 Myr). We estimated this fraction by first integrating the observational mass 
function (Fig. B| middle panel) in the (4-13 Mj»), brown dwarf (13-75 Min) and stellar 


(0.075-10 Mo) mass regimes and then computing the ratio between |FFPs| and brown dwarfs plus 


stars. While previous studies reported a similar fraction of|FFPsf4) our sample doubles the absolute 
number of |FFPs|in a single association, which significantly reduces the statistical uncertainties in 
the mass function. A recent photometric study in the central 6 deg” region of USC}! found a ratio 
of 1.0-1.5 planetary-mass members per square degree, this result is broadly consistent within the 
statistical uncertainties with the ratio we find in the 171 deg? area covered by our survey (0.5-1.0 


planetary-mass members per square degree). 


In the high mass regime (> 1 Mo), our mass function has a slope of [ = —1.2 + 0.2 
(measured in units of logarithmic mass), compatible with the Salpeter slope2 and with all the 
models and simulations considered here 48155154 In the substellar mass regime (< 75 Mj), our 
observational mass function (Fig. [3| bottom) has a slope of T = 0.62+9:}3 (measured in units of 
logarithmic mass), similar to values reported in the field population for L, T and Y dwarfs?314 and 
to models including several channels of substellar object formation£4 The mass function over 
the low-mass stars and high-mass brown dwarfs regime (0.03-1 Mo) is compatible with a log- 
normal distribution! When integrating the analytical mass function of models including mostly 
core-collapse formation Ml over the planetary mass range (4-13 M jup), we find that they predict a 
fraction of only 0.009 — 0.019 underestimating up to seven times our measurement (depend- 
ing on the age assumed). This excess of [FFPs| with respect to a log-normal mass distribution is in 
good agreement with the results reported in o Orionis? Interestingly, our observational mass func- 
tion also has an excess of low-mass brown dwarfs and |FFPs| with respect to simulations including 
both core-collapse and disc fragmentation®. This suggests that some of the [FFPs] in our sample 


could have formed via fast core-accretion in discs rather than disc fragmentation. We also note 


that the continuity of the shape of the mass function at the brown dwarf/planetary mass transition 


suggests a continuity in the formation mechanisms at work for these two classes of objects. 


Hereafter, we use the current knowledge of planet and star formation to discuss the origin of 


The fraction of observed (ferp observed) is the sum of formed by ejection from a 


disc (frrp ejected), the formed by core-collapse (frp core collapse) and the formed by other 


mechanisms (ferP other, including photo-evaporation and ejection from a prestellar cluster). The 
fraction of ejected from a planetary system (frp ejectea) depends on the fraction of stars and 
brown dwarfs that form giant planets (feiant), on the fraction of such planetary systems that become 
unstable (funstabie), and on the number of ejected planets per unstable system (Nejectea). This for- 
mulation is similar to that of previous studies and has the following expression in our study, for 


objects more massive than 4 M jup. 


Free observed — Free ejected F FFP core collapse ch Free other 


13 MJup (1) 


= Í. giant * Í unstable * Mejected + f Ecore collapse (m) dm J f FFP other 
4 MJup 


The fraction of stars that forms giant planets (foiant) is constrained by the observed demo- 
graphics of giant exoplanets measured by radial velocity, direct-imaging, transit and micro-lensing 
surveys DIS), The planetary occurrence rate depends on many astrophysical parameters (the 
host star and planet masses, orbital separation, stellar metallicity, and others), and each of the 
techniques above mentioned is sensitive to a specific region of the parameter space® Therefore, 


we combined the occurrence rates obtained with different techniques to minimise the possible ob- 


servational biases. In Table |1| we summarise the occurrence rates obtained by different authors 
and describe the properties of each study. The fraction of planetary systems that become unstable 
(funstable) Must be at least 75% to match the observed distribution of giant planet eccentricities4:68l 
A minimum of two planets per system is needed for instability to happen, and indeed, many gi- 
ant exoplanets are found in multiple planetary systems or contain hints (such as radial velocity 
trends) of additional companionsél The number of ejected planets per unstable system (Nejectea) 
scales with the number of planets involved in the instability We consider a simplified scenario 
in which planetary systems contain two to four giant planets, the three cases being equally likely, 
and every time a system becomes unstable it ejects one planet. With these assumptions, we find 


Nejectea ~ 1/3 -2 -1/2 + 1/3 -3-1/3 +1/3-4-1/4~ 1. 


Combining the upper and lower limits on the three factors defining the fraction of 
ejected from planetary systems, we obtain frrp ejectea ~ 0.005 — 0.021, which represents between 
10-130% of the [FFP] population we found. The large uncertainty in this percentage is due to the 
lack of precise masses in our observations on the one hand (which itself is related to the lack of 
precise individual ages) and to the uncertainty on the occurrence rate of planets and their ejection 
process on the other hand. A percentage of planets formed by ejection above 100% is obviously of 
no physical meaning and only reflects the limitations of our assumptions and simplifications, which 
are based on the best current knowledge of planetary systems and the latest evolutionary models. 
This result nevertheless suggests that ejection from planetary systems is a significant mechanism 


for formation since at least 10% of them must have formed by ejection from a disc. 


The discovery of a large population of in and also bears important implica- 


tions on the formation and early evolution of planetary systems and, specifically, on the timescale 
of the processes involved. N-body simulations indeed predict that most dynamical instabilities 
happen within 0.1-1 Myr of the planets” formation, although there do exist configurations that 
produce later instabilitie #6] Our results suggest that giant planet systems must form and be- 
come dynamically unstable within the observed lifetime of the region of 3— 10 Myr to contribute to 
the population of While the instability among the Solar System’s giant planets was much 
less violent than those for the mass range of [FFPs| in our sample! current studies suggest that it 


may have also happened early 212 


Instabilities can also be the result of close stellar encounters: numerical simulations have 
shown that dynamical interactions with other stars in clustered environments may induce instability 
in planetary systems or even liberate planets, enriching the [FFP| population 2174) Recent studies 
showing that the demographics of exoplanets depends on the stellar environment confirm that 
such interactions must indeed play a role and contribute to the observed population of Our 
observations suggest that these encounters might take place within the first 10 Myr of a system’s 


life. 


The combined contributions of from core-collapse (13-118%) and ejection from plan- 
etary systems (10-130%) derived from our analysis can explain the formation of the majority of 
But other mechanisms are known to be at work: photo-erosion of prestellar cored=214 has 


been observed around massive B star I] Since is an OB association, some of the 
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might have formed by photo-erosion. At the same time, given that photo-erosion can only occur 
in the direct vicinity of relatively scarce OB stars, we can reasonably assume that the contribu- 
tion of this mechanism to the overall population must be fairly small. The contribution of 
dynamical ejections of sub-stellar embryos in the parent prestellar cluste2B1 could in principle 
be significant source of However, recent hydrodynamical simulations show that while such 
dynamical ejections can produce a realistic population brown dwarfs, they under-produce 
(see Figure B). This new sample of 1s by far the largest and most comprehensive known 
to date and brings exciting opportunities to better understand their origin by performing statisti- 
cally robust studies of their properties. The multiplicity, kinematics and properties of discs among 
are expected to depend on the processes at work and should hold important clues on their 
formation and early evolution. Finally, this new sample also provides an opportunity to study the 


atmospheres of planetary-mass objects in the absence of a blinding host star. 


Methods 
DANCe catalogue 


The COSMIC-DANCe projec? aims at performing deep and complete censuses of young 
nearby associations down to the planetary mass regime. While the identification of young stars and 
even brown dwarfs has become trivial with the advent of Gaia, finding the extremely faint free- 
floating planets well beyond Gaia’s sensitivity limit remains a difficult challenge that requires deep 
ground-based observations. The strategy chosen for COSMIC-DANCe relies on the measurement 


of proper motions and multi-wavelength photometry and the use of modern data mining techniques 
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to identify the faintest members. We therefore combine optical and infrared wide-field images 
found in public archives with our own observations to obtain multi-epoch and multi-wavelength 
catalogues of young clusters and star forming regions in the solar neighbourhood. The relatively 
long time base-line encompassed by this data set (~ 20 years) allows us to derive proper motions 


with a precision of the order of 1 mas yr”!. 


We searched in the European Souther Observatory (ESO) 
and 
public archives for wide-field images inside 


the area 


235° < RA < 252°, 
(2) 
—29.5° Dee < —16.7°. 


We complemented the data found in these public archives with our observations with the 


mounted on the Blanco telescope at the 
the VISTA and VST telescopes at the MegaCam camera at 
France-Hawaii Telescope (CFHT)| the NEWFIRM camera mounted on the 4 m telescope at the 
IRC Peak National Observatory (KPNO}] and (CTTO) the [Hyper Sapam Cam (ASC) mounted on 
the Subaru telescope and the [Wide Field Camera (WFC)| mounted on the [Isaac Newton Telescope] 


(INT), Supplementary Table 1 gives an overview of the various instruments used for this study. 


We used the MAXIMASK and MAXITRACK softwares”!to detect problematic pixels (cosmic rays, 
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dead/hot pixels, satellite trails, saturated stars and associated blooming artefacts, diffraction spikes) 
and problematic astronomical images (e.g. tracking lost). Additionally, we discarded a few images 
after a visual inspection because of their poor quality, limited sensitivity, or acquisition problems. 
Finally, we collected 80 818 individual images of 18 different instruments, obtained over the past 
20 years. The raw and processed data added up to almost 120 TB and were processed and analysed 
on a dedicated HPC server. We extracted over 1.3 billion individual detections from these images. 
The image reduction and the photometric and astrometric analysis are described in a previous 
article=2 Briefly, all images were either processed using their official pipeline (e.g. 
or using Alambic® following standard procedures. Sources were detected, extracted and 
their astrometry and photometry measured using SEXTRACTOR®! and PSFEX2 The astrometric 
and photometric calibration was then obtained using SCAM! and nightly sets of individual 
images were average-combined (weighted by their exposure times) using SWarpl Sources were 
then extracted in these deeper nightly stacks, and proper motions and photometry measured using 
SCAMP again. Deep stacks were also produced combining all available images obtained in a 
given camera+filter and used for the photometry only. The details of the whole procedure can 
be found in the original COSMIC-DANCE article In Supplementary Figure 1 we show the 
precision of proper motions as a function of magnitude. Because of Gaia superiority compared to 
our ground based observations we always use the astrometry from [Gaia DR when available and 


the astrometry from DANCe elsewhere. This explains the shift in precision observed ati ~ 21 mag. 


Our catalogue contains proper motions and photometry (grizyJ H K) for 40 882 164 unique 


sources. To optimise the number of sources with complete photometry (essential for the member- 
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ship analysis, see description of the membership analysis) and accelerate the computational time 
of the membership algorithm, we selected the area where the coverage of most instruments was 
best. We defined an elliptical area centred in (RA = 243.5”, Dec = —23.1°) with a semi-major axis 
of 8.5° in RA and a semi-minor axis of 6.4° in Dec. This selection roughly follows the coverage of 


the UKIDSS near-infrared survey which we also used to define our and surveys. 


The globular cluster NGC 6121 (y, = —12.48 mas yr“, us = —18.9 mas yr”? 6%) is in- 
side the area covered by our survey. To avoid contamination from its members in our sample, we 
discarded the sources encompassed in a circular region of 12’ around the globular cluster centre 
(RA= 245.896°, Dec= —26.527°). The final catalogue contains 28062542 sources and has a 
median precision of < 1 mas yr”? in proper motions for sources brighter than ¿ < 20 mag. Supple- 
mentary Figure 2 shows the density distribution of sources as a function of magnitude for different 
filters. We used the maximum of this distribution as the approximate completeness limit in each 


band. However, it depends on dust extinction and varies with position in the DANCe catalogue. 
Membership analysis 


We used a maximum likelihood approach to infer the parameters of the models describing 
both the cluster and field populations=!!87| The algorithm models first the distribution in the space 
of observables (parallaxes, proper motions and photometry) of the sources that belong to the field 
and then iteratively searches for a maximum likelihood solution for the parameters of the model 
that describes the distribution of observables for the cluster sources. In each iteration, the algo- 


rithm calculates the membership probabilities using Bayes’ theorem and the fractions of cluster 
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and field members as priors. The sources with missing data cannot be used to infer the cluster 
and field models and, for that, 1t is of uttermost importance to define an adequate representation 
space, i.e. a set of parameters which is the largest but at the same time contains a large fraction 
of sources with complete observations. We note however that our algorithm uses the final model 
(computed with complete sources) to obtain a membership probability for the sources with partial 
information by marginalising over the missing information. We searched for members in three dif- 
ferent catalogues: the DANCe catalogue (produced in this work), the Gaia DR2 catalogue and the 
Hipparcos catalogue. These catalogues include very different photometric bands, and we decided 
to run an independent analysis for each catalogue. The parameters (photometry and astrometry) 


used in each case are described in the following paragraphs. 
Initial members 


We compiled a list of 2865 published candidate members in the literature 243188] in the 
area covered by this study. We cross-matched this list with each of our three catalogues (Hipparcos, 
Gaia and DANCe) to obtain the initial list to start each analysis. In the case of Hipparcos, we 
excluded Antares (o Sco) since it is a giant star, and therefore, it falls out of the empirical pre-main- 
sequence isochrone. For the analysis with Gaia and DANCe, we excluded the most extinguished 


members since they confuse our empiric isochrone. 
Representation space 


For the analysis with Hipparcos, we searched for members in the space of pmra, pmdec, 


parallax, V, B—V, where all the sources in the catalogue have complete observations. For the 
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analysis with Gaia, we used the representation space pmra, pmdec, parallax, Grp, G — Grp» 
excluding the Cor band which is less accurate for cool dwarfs2 In this space, 7768 856 sources 
(97%) have complete observations. For the analysis with the DANCe catalogue, we used the rep- 
resentation space pmra, pmdec, 1, J, H, i — Ks. We combined the ¿ band in the optical, which 
has the largest coverage, with the infrared bands J, H, Ks where the ultracool dwarfs are best de- 
tected. With this representation space, 10 483 667 sources have observations in all the photometric 
bands, which represents 37% of the catalogue. We decided not to include the g, r, z, and Y bands 
in the representation space because they reduced the number of sources with complete photometry, 


specially for the coolest objects. 


Field model 


The model of the field population is a[Gaussian mixture model (GMM)|in the whole repre- 


sentation space. We explored models with different number of Gaussians and used the 


information criterion (BIC)|criterion to chose the final model. Since the Hipparcos catalogue has 


a reduced number of sources, we explored models with several components between 1 and 20 and 
selected six as the optimum choice according to the [BIC| For the Gaia catalogue, we explored 
models between 20 and 180 Gaussians and chose 60 as the optimum choice according to the [BIC] 
Finally, for the DANCe catalogue, we explored models between 60 and 300 and chose 100 com- 


ponents. 
Cluster model 
The proper motion distribution of the region of and is much more complex than 
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that of open clusters previously analysed with this methodology 8790] While open clusters have 
symmetric nearly Gaussian distributions in astrometry, this young region shows a rich substructure 
far from Gaussian (see Supplementary Figure 3) and indicative of multiple kinematic populations. 
To model this complex distribution, we used a where the Gaussians are not necessarily 
concentric and explored models with between 1 and 10 Gaussians. Since the Hipparcos catalogue 
contains a very reduced number of sources, we found that a single multivariate Gaussian function 
suffices to model the cluster proper motions and parallaxes. The Gaia and DANCe catalogues are 
much larger and the number of Gaussians selected according to the [BIC|criterion is 5 — 7 for Gaia 


and 5 — 6 for DANCe (depending on the parameter pin, see below). 


We ran the model with different internal probability thresholds®4 i.e. different degrees of 
completeness and contamination (pin = 0.5, 0.6, 0.7, 0.8, and 0.9), and for each we computed the 
optimum threshold, p,p,, using synthetic datt UE, In Supplementary Table 2, we show Pin, Popt> 


and the number of members for each independent analysis (Hipparcos,|Gaia DR and DANCe). 
Final list of members 


Membership probabilities obtained from the analysis with different p;,, values, as well as the 
astrometry and photometry used in the Hipparcos, Gaia, and DANCe catalogues, are available at 
Choosing the best solution (the best pin) is 
a non-trivial decision, which depends on the aim of the study. To study the magnitude distribution 
and mass function, we need a list of members as complete as possible. For this reason, we prefer 


solutions with low Pin values which have greater completeness, although they can also be slightly 
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more contaminated. 


First, we compared the Gaia solutions obtained with different pin and found that 2603 
sources (94%) are the same in all the lists. Additionally, the contamination rate and true positive 
rate computed with synthetic data are very similar in the five studies (see Supplementary Table 2), 
so we had no prior reason to prefer one list to another. Therefore, we chose the list of p;, = 0.5 
as the final list of Gaia since it was the one with the largest number of members. Following an 
analogous procedure with the DANCe solution, we also chose the list of pin = 0.5. With the Hip- 
parcos study, we selected the solution of Pin = 0.7, which represents a good compromise between 
low contamination and high completeness. To this final list, we added the giant star Antares manu- 
ally. Our final list of members contains 3 455 sources from the Hipparcos, and DANCe 


catalogues. 
Membership completeness 


The completeness of our membership analysis depends on the completeness of the astro- 
photometric catalogue and the membership algorithm. Our optical and near-infrared 
VISTA images and archival UKIRT images cover the entire area and ensure that the instrumental 
1, Y, J, H, Ks sensitivities are fairly homogeneous spatially. To get an estimate of the complete- 
ness in the substellar mass regime, we propagated the apparent magnitude completeness of the 
DANCe catalogue to masses. The limiting magnitude to search for ultra-cool dwarfs is set by the 
1 band, which we estimated to be sensitive up to ¿ ~ 26 mag and complete up to 1 ~ 23 mag (see 


Supplementary Figure 2). This approximate magnitude limit completeness corresponds to masses 
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between 7 M jup (assuming an age of 3 Myr) and 13 M jup (assuming an age of 10 Myr). However, 
our membership algorithm is expected to miss highly extincted objects (Ay = 3 mag) and sources 


with near-infrared excess related to the presence of circumstellar material. 
Membership validation 


The membership classification mostly coincides in the Hipparcos—Gaia and Gaia-DANCe 
studies, in the magnitude range where both catalogues are complete. The small differences be- 
tween catalogues can be attributed to the different information provided by each of them. We 
used the Gaia membership analysis, with the additional information on the parallax, to estimate a 
contamination rate of 8% on the DANCe membership in the magnitude range 9 < J < 14 mag as- 
suming that the Gaia selection is perfectly clean. Similarly, the comparison between the Gaia and 
DANCe samples over the common luminosity domain shows that one third of the objects identified 
with Gaia are not recovered with DANCe because of either missing photometry, high extinction 
or near-infrared excess likely related to the presence of a circumstellar disc. Therefore, the com- 
pleteness of our census is expected to be better in{USC]than in[Oph] since extinction is much lower 
and near-infrared excesses related to discs should be less frequent given the more advanced age 


and timescale for disc decay. 


We note that at J ~ 10 mag (see Fig. |2) there are some contaminants. These are sources 
identified with the DANCe membership analysis, therefore, using only the proper motions and 
photometry (no parallaxes). The proper motions are compatible with the proper motion distribution 


of the association and that is why they are classified as members. Additionally, since they are 
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bright stars, their photometry saturates in many bands so we have little photometric information. 
In any case, these represent a 1% of the members which is within the contamination rate that we 


estimated. 


We recovered the majority of members previously reported in the literatur 12:45] Supple- 
mentary Figure 4 shows a comparison between the members reported by previous studies and 
the members found in this study, using the photometry and proper motions we measured in this 
work. In this Figure, we only considered studies sensitive to[FFPs| We missed around 80 substel- 
lar members previously reported in the literature, most of which are in [Oph] are highly extincted 
or host circumstellar discs (as can be seen in the colour-magnitude diagrams). Besides, some of 
the members reported in the literature are discarded by our membership analysis because of their 
inconsistent proper motions. This is especially significant in a recent study“ where only half of 
their members are identified in our astro-photometric analysis. The other half are either classified 
as non-members (having photometric and/or proper motion measurements inconsistent with the 


association) or undetected. In this study, we add ~800 of new members, 70-170 of which are 


depending on the age assumed. 


The mass function 
Distances 


We used Kalkayotl|™ to infer Bayesian distances for all the members with a parallax mea- 


surement in the |Gaia DR2|catalogue. We used a Gaussian prior with a locus and scale of 145 pc 


https://github.com/olivares-j/kalkayotl 
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and 45 pc, corresponding to the median and three times the standard deviation of the distribution 
of distances obtained inverting the parallax. For the sources in the DANCe catalogue, without 
a parallax measurement, we sampled the distance from the cluster distance distribution obtained 


with all the Gaia members. 
Masses 


We combined the apparent photometry (grizy JH Ks) and the distance estimate of each star 
to obtain absolute magnitudes. These were compared, in a Bayesian framework, to theoretical evo- 
lutionary models to infer the posterior distribution of the mass and extinction of each source with 
Sakar". This algorithm ignores any possible source of error related to the theoretical evolution- 
ary model chosen by the user. The model does not include effects on the variability of the source 
due to binarity, activity, or other factors. These effects eventually end up included in the extinction 
estimate, enlarging its uncertainties. Finally, the mass and the extinction are degenerated but at 
least half of the planetary-mass objects we found are in regions of low extinction, favouring their 


planetary nature. 


We had to use different models for the high and low-mass regimes as there is no single 
set of models covering the entire mass range of our members. We combined the PARSEC- 
COLIBR]} models and BHACIi models“ which cover the high and low-mass range of our 
members, respectively. We find that both grids of models agree fairly well around 0.5 Mọ and de- 


cided to use the masses inferred from the BHAC15 models below 0.5 Mo and the masses inferred 


https://github.com/olivares-3/Sakam 
*http://stev.oapd.inaf.it/cgi-bin/cmd 


“http://perso.ens-lyon.fr/isabelle.baraffe/ 
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with the PARSEC-COLIBRI models above 0.5 Mo. In Supplementary Figure 5 we show two ex- 
amples of the best-fit spectral energy distribution models obtained with Sakam for a brown dwarf 
and a planetary-mass object. The complete final list of members with the masses and extinctions 


inferred with Sakam is available at the[CDS 
Magnitude and mass distributions 


To obtain the magnitude (mass) distribution, we sampled the individual magnitude (mass) of 
each source with a Gaussian centred at the measured magnitude (mass) and a standard deviation 
equal to the uncertainty. Then, we defined a grid between the least and most bright (massive) ob- 
ject in our sample and added the contribution of all the sources to each magnitude (mass) bin. We 
convoluted this distribution with a Gaussian [kernel density estimation (KDE)] with a bandwidth 
chosen according to Scott’? and Silverman’ rules. We estimated the uncertainties in the mag- 
nitude (mass) function with a bootstrap of 100 repetitions and reported the 1o and 30 confidence 


levels. 


Data availability The data that support the findings of this study will be available at the CDS 


after the reviewing process. 
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Figure 1 Sky distribution of stars (gold triangles), brown dwarfs (blue squares), and 
(red dots) discovered in this study and classified assuming an age of 5 Myr. The 
dashed ellipse indicates the area analysed with the DANCe catalogue (see methods). 
The background images are in the optical (credit: Mario Cogd%) and at 857 GHz (credit: 
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Figure 2 Colour-magnitude diagram of the members of|USC|and|Ophlidentified in this 
work: previously known members (gray) and new members (black). The error bars 
represent the uncertainty in the photometry reported in the Gaia and DANCe catalogues. 
The BHAC15 isochrones“ (solid lines) and the PARSEC-COLIBRI isochrones“! (dashed 
lines) at 3 Myr (red) and 10 Myr (blue) as well as the extinction vector are overplotted. 
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Table 1: Planetary occurrence rate for different studies. 


Method POR (%) Mass Separation | SpT star FFFP ejected Zo, 

R 2.1 +0.5 4-13M; | 0.1-100 AU | solar type | 0.016-0.021 | 23—130% 
DÉI 0.6187 5-13 M; | 30-300 AU | BAFGKM || 0.005-0.006 | 7-38% 

D Led 5-13 M; | 30-300 AU | BAFGKM | 0.014-0.018 | 20-115% 

MLS® 1—2 | 0.007-0.02* | 0.2-5 Op all 0.008-0.02 | 11-125% 

RV+DI+ML95 ~4 1-13 M; | 2-1000 AU all 0.03-0.04 | 44-250% 


Col. 1: planet detection method: radial velocity (RV), direct imaging (DI), microlensing 
(ML), Col. 2: planet occurrence rate (POR), Col. 3: mass range of the planets, Col. 4: 
separation or semimajor axis, Col. 5: spectral types of the primary body, Col. 6: estimated 
fraction of ejected FFPs using Eq.[t]and assuming funstable = 0.75 — 1 and Nejectea = 1 (See 
main text), Col. 7: estimated percentage of observed FFPs which were ejected from 
planetary systems. The conservative ranges we provide include both the uncertainties 


from our observations and current knowledge of planetary systems. 


* mass ratio between the planet and star mass (q). 
t our observations are not sensitive to masses between 1-4 M; and thus, the POR is 
overestimated with respect to our observations. 
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Figure 3 J apparent magnitude distribution (top) and mass function (middle and 
bottom) of the members of[USC]and[Oph| The shaded regions indicate the 1 and 30 
uncertainties from a bootstrap (top) and the dispersion due to the age (3-10 Myr, middle 
and bottom). The mass functions from simulationg/8 are overplotted on our 
observational mass function (bottom). All the functions are normalised in the mass range 
0.004—10 Mo. The hydrogen (75 My,,,) and deuterium (13 M.;,,,) burning limits are 
indicated by the vertical dotted lines according to the BHAC15 evolutionary models and 
assuming an age of 5 Myr. 
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Supplementary Figure 1: Estimated proper motion error as a function of i magnitude for the 
DANCe catalogue. We note that whenever there is a proper motion measure from Gaia we use 
it. This explains the change in precision observed at 7 ~ 21 mag. 
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Supplementary Figure 2: Density of sources per 0.2 magnitude bin as a function of magnitude for 
all the sources in the DANCe catalogue. 
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Supplementary Figure 3: Vector point diagram (left) and parallax—proper motion diagram (right) of 
USC and Oph. The members are shape- and colour-coded according their origin: Gaia and DANCe 
analysis (red circles), only Gaia analysis (blue left-pointing triangle), and only DANCe analysis 
(green right-pointing triangle). The error bars represent the astrometric uncertainties reported in 
the Gaia and DANCe catalogues. 
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Supplementary Figure 4: Vector point diagrams (top) and colour-magnitude diagrams (bottom) 
of the members found in this study (grey dots) and the members reported in previous studies! 
(coloured markers). The members in common are indicated by a black edge line. The error bars 
represent the observational uncertainties in the Gaia and DANCe catalogues for the members pre- 
viously identified in the literature and discarded by our membership analysis. We note that our 
study is not sensitive to highly extinct areas or the detection of sources with moderate/significant 
near-infrared excess. 
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Supplementary Figure 5: Spectral energy distribution of two of our substellar objects: a brown 
dwarf mass object (left) and a planetary-mass object (right). The photometric observations (black) 
and the dereddened observations (blue) are indicated. The error bars represent the photometric 
uncertainties and are smaller than the markers in most cases. The BHAC15 models at 5 Myr used 


for the best-fit models are shown (gray). 
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Supplementary Table 2: Performance of the membership analysis. 


Hipparcos Gala DANCe 

Pin | Dopp Memb CR TPR) Po: Memb CR TPR | pp Memb CR TPR 

4) (Gi Eo) (Gi (%) (Gi 
0.5 (0.97 116 17 92 |0.95 2762 1.1 99.2 | 0.77 2556 2 98 
0.6 | 0.96 108 7 95 |0.96 2698 0.9 99.2 | 0.78 2458 2 99 
0.7 || 0.94 112 7 96 |0.96 2678 1.0 994 [0.83 2342 15 99 
0.8 || 0.95 103 5 96 |0.96 2661 0.9 996/088 2185 09 99 
0.9 || 0.96 78 3 98 |0.96 2623 09 996 |0.83 2086 0.9 99 


Performance of the membership analysis obtained with different internal probability 
thresholds (p;,) for the three catalogues considered, namely Hipparcos, Gaia, and 
DANCe. For each internal probability threshold (pin) we show the corresponding opti- 
mum probability threshold (p+), number of members (Memb), contamination rate (CR) 
and true positive rate (TPR). The pop, CR, and TPR were obtained with synthetic data 
(see text). 


10 


